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Abstract 

The Morris-Lecar (ML) model has applications to neuroscience and cognition. A simple 
network consisting of a pair of synaptically coupled ML neurons can exhibit a wide 
variety of deterministic behaviors including asymmetric amplitude state (AAS), equal 
amplitude state (EAS), and steady state (SS). In addition, in the presence of noise this 
network can exhibit mixed- mode oscillations (MMO), which represent the system being 
stochastically driven between these behaviors. In this paper, we develop a method to 
specifically estimate the parameters representing the coupling strength (<? S yn) and the 
applied current (i a pp) of two reciprocally coupled and biologically similar neurons. This 
method employs conditioning the likelihood on cumulative power and mean voltage. 
Conditioning has the potential to improve the identifiability of the estimation problem. 
Conditioning likelihoods are typically much simpler to model than the explicit joint 
distribution, which several studies have shown to be difficult or impossible to determine 
analytically. We adopt a rejection sampling procedure over a closed defined region 
determined by bifurcation continuation analyses. This rejection sampling procedure is 
easily embedded within the proposal distribution of a Bayesian Markov chain Monte 
Carlo (MCMC) scheme and we evaluate its performance. This is the first report of a 
Bayesian parameter estimation for two reciprocally coupled Morris-Lecar neurons, and 
we find a proposal utilizing rejection sampling reduces parameter estimate bias relative 
to naive sampling. Application to stochastically coupled ML neurons is a future goal. 
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1 Introduction 



Transmembrane voltage is often recorded during physiological study of biological neurons. 
However, voltage-gated ion channel activity and neurotransmitter levels are quite difficult 
to measure directly and are usually unobserved in such studies. In addition, there is a great 
diversity of neuron morphology, protein expression, and plasticity which may affect voltage 
dynamics and synaptic transmission (DeCarli et al., 20f2; Kollins and Davenport, 2005). 
Early development and senescence may also be major determinants of voltage response pro- 
files (Yeoman et al., 2012; Liu et al., 2012). Synaptic tuning in particular is thought to be 
an essential mediator of learning, stimulus response integration, and memory. There is evi- 
dence that memory and learning may depend critically on several distinct types of dynamic 
behavior in the voltage of neurons. 

The ML model reproduces the voltage of a single neuron and, depending on parameterization 
and initial conditions, can exhibit many of the experimentally observed behaviors of biolog- 
ical neurons (Morris and Lecar, 1981). In this paper, we explore a simple neural network 
consisting of two biologically identical, reciprocally coupled ML neurons. Yu et al. (2008) 
have shown that this modest model can exhibit a wide range of oscillating or non-oscillating 
voltage depending on the values of just a few parameters, specifically in this study, / app and 
g syn . In the absence of noise, the model can predict synchronous or asynchronous firing, as 
well as either equal or unequal action potential amplitudes. Additionally, in the presence 
of even small noise in the applied current and weak synaptic coupling, the system can ex- 
hibit mixed-mode oscillations (MMO) characterized by periods of small amplitude oscillation 
interrupted by large amplitude excursions. 

In further work with the two ML neuron model, Thompson (2010) explored two synaptically 
decoupled neurons driven by both common and independent intrinsic noise terms. They 
found that shared common noise promotes synchronous firing of the two neurons, while 
separate intrinsic noise terms promoted asynchronous firing. The relative scaling of the two 
noise sources was observed to be key in predicting the degree of synchrony. In addition, 
while they did not specifically look at MMO, they hypothesized that such synchrony in a 
synaptically coupled network would increase the probability of MMO, by facilitating longer 
residence times within the unstable periodic orbits adjacent to the system's stable periodic 
orbits. Indeed, in this paper we will detail the relative positions of these parameter regions 
as they are of key importance to our conditioned likelihood approach. Specifically, we will 
provide a quick look-up table for the region in parameter space where stable periodic orbits 
are possible. 

Ditlevsen and Samson (2012) develop a expectation-maximization (EM) stochastic particle 
filter method to estimate the parameters in a single ML neuron based on observation of 
voltage only. A key aspect of their approach is that they assume both the voltage and the 
channel gating variables are in an oscillatory regime, but stochastically perturbed. These 
perturbations are considered nuisance parameters which their method marginalizes away. 
Specifically, they treat the unobserved channel gating variable from the model as a com- 
pletely latent variable. Starting from estimates of the initial conditions for the voltage and 
channel gating variables, they iteratively predict the gating variable and voltage and then 
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update the predicted voltage to the next time step using a modification of the well-known 
Euler differential equation solver. They discuss that an assumption of stationarity in their 
method limits applicability to only short time windows over which current input can be 
considered constant (e.g. 600ms). They also note that certain parameters, conductances 
and reversal potentials in particular, are sensitive to choice of tuning parameters required 
by the method. 

These studies demonstrate the active progress as well as the challenges of model parame- 
ter estimation for biological neuronal models and, more generally, for relaxation oscillator 
models. Each of these studies derives asymptotic approximations or general forms for model 
likelihood, but use fundamentally different techniques and assumptions in doing so. In each 
study the approach is specifically crafted to the model. In this paper we attempt to develop 
a convenient Bayesian estimation scheme with only a few tuning parameters and relatively 
few mild assumptions. We focus our attention on deterministic synaptically coupled ML 
neurons. Application of our method to stochastically coupled ML neurons is on-going work 
in our group. 

In the case of ML, estimation of / app and g syn is non-trivial due to the diversity of possible 
dynamic behavior and the abrupt transitions among these seen with just small changes in 
these parameters' values. However, we can better understand the critical values of these 
parameters by studying the system's bifurcation structure. We are able to locate parameter 
regimes where dramatic changes in the system appear. The neurons analyzed in this study 
are classified as Type II neurons, characterized by discontinuous drastic shifting between 
behavioral states. Because there is a distinct switch in behavior, bifurcation analyses de- 
termine a closed region of parameter space over which the relevant dynamics may occur. 
Sampling over such a feasibility region amounts to conditioning the inference on an a priori 
assumed class of dynamics (e.g. stable node, limit cycle, steady state etc.). While facilitat- 
ing conditioning the likelihood on feature statistics of the voltage, this may translate into 
increased confidence and reduced bias in the parameter estimates. 

2 Reciprocally Coupled ML Model 

Our goal is parameter inference based on the temporal voltage response of two synaptically 
coupled neurons which are deterministically coupled to voltage-gated ionic conductance dy- 
namics (Morris and Lecar, 1981). A single ML model has a two-dimensional phase space 
and is known to reproduce many of the behaviors experimentally observed in biological neu- 
rons (Yu et al., 2008). Therefore, systems of coupled ML neurons may offer a reasonable 
starting point for developing statistical inference methods for models of neuronal networks. 
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The ML network we study is, 

-^r = ^(fi'ca • moo (vi) ■ {vi - v Ca ) - Qk ■ wiivi - v K ) - 9l ■ wi{vi - v L ) 

+ 4pp - g S yn ■ Si • (Vl - V syn )) + 5-£i 

= ^(fi'Ca • moo (v 2 ) ■ (v 2 - v Ca ) - g K ■ w 2 {v 2 - v K ) - g L ■ Wi{v 2 - V L ) 

+ 4pp - fi'syn • S 2 ■ (V 2 - V syQ )) + 5 ■ £ 2 

dsi _ Sqq (v 2 ) - si 
dt t 

dS 2 _ Sqq (fj) ~ S 2 

dt t 
dw 1 



where 



dt 
dw 2 

~dT 



= A(ul)(woo (vl)-wl) 
= A (v2) (w M (v2) - w2) 



If fx — V 11 

m - (X) = 2 V + tanh [-02- 



x—vt 

1 + e vs 



. . 1 / /rr — f 3 

w„o (rr = - 1 + tanh 



2 V \ v4 

x — v3 



A (x) — 6 cosh 



Note that in the stochastic version of this model 5 > and £1 and ^ 2 are standard independent 
Wiener process variables. In this paper, however we will be concerned with the deterministic 
version of this model where 5 = 0. 



Table 1. Parameters of ML model 

Assumed values of these parameters are given in the Appendix. Initial conditions of the variables may also 
determine the observed dynamics. 



Variable 


Definition 


Variable 


Definition 


-fapp 


Applied current 


^1,^2,^3,^4 


Membrane potentials 


9syn 


Coupling strength 


i>L, »Ca> VK, fsyn 


Equilibrium potentials 


c 


Membrane capacitance 


mirij> 


(1 + tanh[(y - V 1 )/V 2 ])/2 




Conductance of membrane channels 


Sl, S2 


|.(l + tanh[^]) 


W1,W2 


Recovery variables 





While the many parameters of the ML model are, in principle, experimentally verifiable, they 
impart a high dimension parameter space determining the behavior of the system. Indeed, 
there is a non-trivially large diversity of behaviors already possible by varying just a small 
few of these parameters. To simplify our exposition and make the essence of our approach 
clear, the only two parameters we examine in detail are J ap p, the exogenously applied current, 
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and g syn , the synaptic coupling strength. Hence, in this study, 7 app and g syQ are the only 
unknown parameters and other parameters are assumed to have the values given in the 
Appendix. In order to estimate these parameters, we develop a Bayesian MCMC method 
based on a Metropolis-Hastings sampling approach with conditioned likelihood and rejection 
sampling proposal distribution. 



3 Bifurcation Continuation Analyses 

Often a slight change in the parameters of a system causes a characteristic alteration in the 
system's behavior. The study of this change in behavior is known as bifurcation analysis. 
Information about the behavior of the ML model can be obtained by studying bifurcation 
diagrams. In particular, parameter values corresponding to steady state and periodic solu- 
tions can be found. This information can then be used to better approximate J app and g syn 
using MCMC. 

For a single neuron, the bifurcation from steady state to oscillation is diagnostic for two broad 
classes referred to as Type I neurons and Type II neurons. In Type I neurons (Figs, la) and 
lb)), the switch from a SS to an oscillatory state is gradual and can often be difficult to detect 
precisely. However, in Type II neurons (Figs, lc) and Id)), the change in state is sudden and 
drastic with a discontinuous jump in the frequency of action potentials. The neurons used 
in our study are determined to be Type II neurons. The sudden switch in behavior allows 
us to refine a proposal distribution for MCMC based on candidate parameter values. 
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Figure 1. Type I vs Type II neurons 

Steady states are black, with a solid line depicting a stable state and a dashed line representing an unstable 
state. Oscillatory states are either green or blue. Green indicates a stable periodic state while blue represents an 
unstable periodic state. All were generated with a g sy „ value of to simulate uncoupled behavior. In a)-b), a 
Type I ML neuron is illustrated where all parameters are same as in Appendix except 113 = 15 and V4 = 15. As 
seen in b), increasing I 3pp continuously increases the firing frequency from zero. In c)— d), a Type II ML neuron 
is illustrated where all parameters are identical to those in Appendix. In contrast to a Type I neuron, there is a 
discontinuous jump in frequency seen in d). 



When two ML neurons are coupled, there are in addition two different types of periodic 
behavior. The system is in an asymmetric amplitude state (AAS) when the two voltages are 
oscillating with different amplitudes. In particular, one neuron experiences large amplitude 
oscillation (LAO) while the other experiences small amplitude oscillation (SAO). Alterna- 
tively, an equal amplitude state (EAS) occurs when both neurons are oscillating with the 
same amplitude of voltage. Examples of all these states can be found in Fig. 2. 
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Figure 2. Stable dynamics of coupled deterministic ML model 

Voltage of neuron 1 is shown in blue and of neuron 2 is shown in red. In (a), SS dynamics develop after an 
initial transient. Parameters are same as in Appendix except that I 3pp = 95.5 and g syn = 0.15. Anti-phasic AAS 
arises in (b), such that neuron 1 is experiencing SAO while neuron 2 is experiencing LAO. Parameters are same 
as in Appendix except that / app = 97.5 and g syn = 0.15 and the initial conditions {«i,«2,M)i,ro2,si,S2} ={- 
3,-20,0,0.17,0,0}. EAS emerges in (c) where voltages of neuron 1 and neuron 2 fire with equal amplitude. 
Parameters and initial conditions are same as for b) except g syn = 0.2. Here the voltages gradually become 
entrained in-phase. 

If the dynamics of a given data set can be characterized (e.g. AAS, EAS, or SS), then a 
corresponding parameter regime can be embedded within the MCMC proposal distribution 
and implemented by simple rejection sampling. In this study, we refine the parameter regime 
purely based on whether the system is in an oscillatory state, regardless of whether it is AAS 
or EAS. Bifurcation analyses were performed computationally using the AUTO library via 
the software XPPAUT (Ermentrout, 2002). Evaluations of the model were performed by 
XPPAUT as needed directly from within MATLAB via the functionality provided by the 
xppauttools package of the Snifnib JAVA library sourcef orge .net/projects/ sniff lib. All 
other aspects of analysis were implemented in MATLAB. 

4 MCMC Estimation Method 
4.1 Metropolis-Hastings Sampling 

MCMC is an iterative random walk method that avoids the need for a closed form of the pos- 
terior distribution function. Rather than a direct method which would require evaluation of 
a high dimensional integral of the posterior over the parameter space, MCMC behaves much 
differently, instead returning a sample of parameters from the posterior. Since these param- 
eter samples are drawn from the desired posterior, consistent estimates of moments such as 
means and variances of the parameters may be easily calculated. The Metropolis-Hastings 
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sampler was used which in lieu of an exact posterior substitutes a proposal distribution. We 
next discuss the proposal, likelihood, and prior distributions supplied code to the MCMC 
routine. 

The first of these is the proposal, which in Metropolis-Hastings sampling supplies new can- 
didate parameter sets based on the current set in the MCMC chain. Many of the ML 
parameters are constrained to be non-negative and in particular we assume this for g syn 
and iapp. However, we accomplish unconstrained parameter estimation over this support by 
way of log transformation. To start, the program is given a set of initial guesses for the 
log transformed parameters I app and g syn . This set of parameters is labeled as 9 . A new 
candidate parameter set is drawn from a Gaussian distribution with mean centered on these 
log transformed parameters, 9*. The standard deviation is a fixed and predetermined value 
typically referred to as a "mixing value." Then, if the quantity resulting from multiplying 
the model likelihood by the model prior for the inverse transformed is greater than that 
obtained for Q , then the candidate parameter set is "accepted" as the new 9 . Otherwise, 
the new set 9*, is "rejected" and 9 remains the initial guess. A new parameter set candidate 
is drawn at random from the same distribution, continuing this cycle until a 9* is "accepted." 
The method is continued until after a burn-in period and the chain is determined to have 
settled into an steady state equilibrium. 

4.2 Construction of Conditioned Likelihood 

Now we discuss implementation of the likelihood. Analytic and/or efficient forms for the 
general voltage likelihood for coupled ML are intractable and unavailable. Instead, we con- 
struct conditional likelihoods based on feature statistics (Fraser, 1964, 2004; Ghosh et al., 
2010). The first feature statistic we consider is the cumulative power. 

4.2.1 Cumulative Power 

For a twice differentiable signal m(t), cumulative power P(t) is defined as 



for real finite t. P(t) is a useful feature statistic for a variety of models. For example, 
square-integrable functions, including bump functions, maxima/minima curves, and satu- 
ration curves all have P(t) G 0(1). In contrast, any finite sum of sines and cosines have 
P(t) G 0(t) (Quinn, 2011). By way of Fourier representation, this 0(t) behavior character- 
izes a broad array of continuous periodic functions (see proof in Appendix). Surprisingly, it 
has even been shown that P(t) G 0(t) in a class of linear stochastic dynamical systems lack- 
ing differentiability at countably many times (Bates et al., 2012). In all such cases, having 
P(t) G 0(t) imparts the intuitive notion that the dynamical system exhibits a consistent 
(stationary) duty cycle and accumulates power (on average) at a constant rate. 
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For periodic oscillating functions such as the voltage in the ML model, the graph of the 
cumulative power is similar to the red or blue line in Fig. 3. 





200 400 600 800 1000 



200 400 600 800 1000 



Figure 3. Cumulative power of voltage 

The blue line is the estimated cumulative power of the observed voltage (Pdata(i))- The red line is the cumulative 
power of the voltage predicted from the ML model (PmodelM) based on candidate parameter estimates sampled 
from the proposal distribution. Dashed lines give the ±1 standard deviation wedge interval. See Eqn.3 for details. 



In Fig. 3, the blue line is the estimated cumulative power of the data set Pdata(^)- The red line 
is the cumulative power of the model P mo de\{t) based on parameter estimates with the dashed 
lines giving the ±1 standard deviation interval. (see Eqn.3 for details). Intuitively, matching 
-fmodei(^) to -fdata(^) imparts qualitative matching in terms of amplitude and frequency and 
may reduce bias, especially in those parameters of the model which determine these features. 



l app 



and g syn strongly affect both of these features. 



The curves were determined by locally weighted polynomial regression (LWPR) which is a 
time-domain smoothing method competitive with frequency domain methods such as But- 
terworth filtering Watkin (2011) and wavelets (Fan and Gijbels, 1996) and are especially 
convenient when the prediction of time derivatives is desired. In standard least squares, the 
over-determined system 



p u px 1 



is taken to have the solution 



bpxi — {X nxp X nxp ) X nxp Y nx i. 



In contrast, in weighted least squares, each equation is weighted by a diagonal matrix 

W nxn 

M^nxn^nxl = ^V nxn X nX pbp X i, 

which is taken to have the solution 



^pxl — (X nX pW nxn X nx ,pj X nX pW nxn Y nx i. 
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Figure 4. Example LWPR on noisy unevenly spaced data 

Noisy data (blue) yi = sin(xi) + £j where £; ~ Gaussian(0, 0.1). Locally quadratic regression using a tricube 
weighting kernel and a nearest-neighbor bandwidth of 53% (red) goes smoothly through the data coordinates. 
The values of the tricube weighting kernel (green) give most weight to the points near the point to be predicted 
(weights for x=5.95 being shown). The locally quadratic regression (magenta) may easily be differentiated to 
generate robust 1st and 2nd derivative estimates at that point. 



The jth diagonal element of W nxn is the evaluation of a weighting kernel having a compact 
support taken at the distance (in time) between the point being predicted and the jth 
datum. Weighting kernels are parameterized by a bandwidth or tuning parameter h which 
determines the distance beyond which the weights go to zero. Fan and Gijbels (1996) gives 
a comprehensive review of typical weighting kernels used in LWPR. In our application, we 
choose a tricube weighting kernel and a nearest-neighbor bandwidth which is illustrated in 
Fig. 4. 



To determine P(t) from LWPR for each datum j G {1, . 
we calculate, 



, n} and for each voltage k G {1, 2} 



Pdata,,k{tj) = 2^ ^k,i f^dlta,fc(^fc,i l^data, 
i=0 L 

Pmodcl,k(tj) = pleak fc Atk t i (^rnodel,fc(^M l^model. 

8=0 L 



/>pow 
L 'k,j 



Gaussian (P da ta(^) - -Pmodei(tj) 0,p0 fc + pscale fc P da ta(ti 



(2) 
(3) 



where Gaussian(x a) is a Gaussian PDF with mean fi and standard deviation a evaluated 
at x. The second derivative estimates are easily obtained via the LWPR (see Fan and 
Gijbels (1996) for details) on each of the voltages. LWPR bandwidths /idata and h mo d e i were 
determined by generalized cross-validation and fixed during MCMC. In the current study, we 
did not estimate but rather fixed p0 fc at the root mean square error from a linear regression 
of -Pdata(^j) on t which gave satisfactory results. By similar reasoning, pscale fc = 1 is used as 
a reasonable initial guess but was estimated during MCMC. Note that in this formulation, 
there is no requirement that the voltages be observed at the same time nor that the spacings 
between times are equal. 
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4.2.2 Mean of the Voltages 



In addition to conditioning on cumulative power, we conditioned the likelihood on the mean 
voltages. Specifically, we calculate, 



/'mean 



Gaussian 



data,fc 



model, k 



mstd 



(4) 



where Vdata is the mean of the observed voltage data for the fcth neuron, Vmodei is the mean of 
the predicted voltage data for the fcth neuron, and mstdfc is the estimated standard deviation 
of the mean voltage residual for the kth. neuron. A reasonable initial guess for the mstdfc is 
lOmV. The conditioned model likelihood is thus obtained by multiplication of the C p k ow and 
C~ likelihoods, 



'mean 
-k 



C{6*) = ul =1 (n? =1 ££7) • ci 

where 9* = j J ap p, g^ n , pscale 1? pscale 2 , pleak 1 , pleak 2 , mstd 2 , mstd 2 |. 



(5) 



4.3 Construction of Parameter Rejection Region 

Bifurcation analysis can be used to learn about the behavior of the system, specifically where 
the system is EAS/AAS or SS. Assuming we know the state of the data set, we can use the 
bifurcation analysis to limit the parameter combinations over which MCMC tests. This 
should lead to a more accurate approximation of J app and g syQ . 

One reason we choose to construct parameter boundaries based on the state of the system is 
because of the vast difference in cumulative power for EAS/AAS versus SS solutions. Since 
MCMC calculates probability based on cumulative power, proposed points in a different 
state than the real data set will always lead to a rejection. In addition, it is easy to observe 
whether a data set is oscillating without knowing the parameter values. 



4.3.1 Limit on g syn 

A system will always exhibit the behavior of the stable solution for a set of parameter val- 
ues. For instance, in Fig. 5(a) for an J app value of 150, there is a stable periodic solution, 
an unstable periodic solution, and an unstable steady state solution present. Depending 
on the initial conditions used, the system may move briefly toward the unstable solutions, 
but eventually it will settle to exhibit behavior consistent with the stable periodic solution. 
Therefore, long-term oscillating solutions will only be present for parameter values corre- 
sponding with a stable periodic branch in the bifurcation diagram. Starting at a g syn of 0, 
it can be observed that the stable periodic branch of the system becomes shorter as g syn 
increases. This can be seen in Fig. 5. The length of the branch decreases between g syn values 
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of 4 and 7 before disappearing altogether at a g syn value of 10. Narrowing this further, we 
find that the value of g syn for which the stable periodic branch ceases to exist is about 9.08. 
At this point, the system is always SS. Therefore, based on our observed behavior, oscillatory 
solutions only occur for g syn values less than 9.08. To maintain biological realism, we also 
impose g syn > which provides a lower bound. 




100 200 300 100 200 300 

lapp(^A/cm 2 ) lapp(nA/cm 2 ) 




lapp(^A/cm 2 ) lapp(nA/cm 2 ) 
Figure 5. Loss of stable periodic orbits 

In a) a there is a wide range of I app values for which stable periodic orbits are possible. In b) the range of I app 
values leading to stable periodic orbits is considerably smaller. In c) stable periodic orbits are not possible and 
(as seen in d) do not return for larger values of iapp. 



4.3.2 Limit on J app 

With an established range for g syn , we wish to constrain the corresponding J app values for 
EAS/AAS oscillations. To do this, we carried out two-parameter bifurcation continuation 
of the limit point (LP) bifurcations. This revealed an envelope enclosing a similar two- 
parameter Hopf continuation. Hopf bifurcations can be of either the sub-critical or super- 
critical type with stable periodic orbits possibly emerging below or above the Hopf bifur- 
cation. We observe the LP continuation envelope extending below the Hopf envelope for 
smaller values of J app and extends above the Hopf envelope for larger values of 7 app . Outside 
of the LP region, SS is the only stable behavior. Table 2 shows the final oscillation parameter 
ranges found. It can be observed that as g syn increases, the LP envelopes draw closer until 
eventually intersecting around g syn =9.08 and / app =95.840. This is consistent with our g syn 
boundary. 
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gsyn 


lapp 


gsyn 


lapp 


0.007 


238.382 


9.019 


95 


842 


0.832 


238.097 


8.199 


95 


867 


1.657 


237.885 


7.380 


95 


889 


2.482 


236.004 


6.561 


95 


909 


3.306 


231.454 


5.742 


95 


925 


4.131 


223.402 


4.923 


95 


938 


4.956 


211.055 


4.103 


95 


945 


5.781 


193.919 


3.284 


95 


946 


6.606 


172.617 


2.465 


95 


939 


7.430 


148.735 


1.646 


95 


918 


8.255 


123.144 


0.826 


95 


871 


9.080 


95.840 


0.000 


95 


724 



Table 2. Proposal rejection look-up table 

These 24 values give good approximation to the region where EAS/AAS oscillations are possible (see Fig. 7). 
Defining a closed region, these few values allow computationally efficient rejection sampling within the proposal 
distribution function. 




Figure 6. 3D bifurcation diagram 

Hopf and LP enclosures (see 2D projection in Fig. 7) are shown intersected by planar bifurcation diagrams taken 
at cross-sections g syn = 4, 7, and 10. 
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4.3.3 Rejection Sa 
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Figure 7. Rejection region for oscillating solutions (AAS and EAS) 

As in Fig. 5a), the Hopf bifurcation continuations (cyan, red, and magenta dash) exhibit sub-critical (red) or 
super-critical (magenta) stable periodic orbits. Accordingly, stable periodic orbits are found for any parameter 
combinations bounded within a slightly larger region (gray) having an area 933.75 ■ mS ■ cm~ 4 and 
coinciding with LP bifurcation continuations (green and blue solid). Outside of this region, SS is the only stable 
behavior. 



Determining if a parameter candidate is within the 2D region specified in Table 2 can be easily 
accomplished by a variety of methods. However, we have interest to extending our methods 
to more parameters than just / app and g Byn . Therefore Delaunay triangulation was adopted as 
a general approach to the in/out region testing. In this approach the possibly D-dimensional 
region is decomposed into simplices (triangles in 2D) having D + 1 vertices such that the D 
dimensional circumspheres about any vertex do not contain any other points (Cignoni et al., 
1998). The number of simplices is determined by the number of points defining the region 
boundaries which is our case result from parameter bifurcation continuation via XPPAUT. 
This triangulation process is efficient and may even be performed offline prior to MCMC and 
thus represents negligible computational expense. As seen in Fig. 8, a parameter candidate 
generated by a random draw from the proposal distribution is determined to be within this 
region if it is internal to any one of this simplices. For a region determined by m continuation 
points, this is an 0(m) operation and quite efficient. 
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lapp (|a,A/cm 2 ) 

Figure 8. Delaunay triangulation of "NO ROLLBACK" region 

In an approach which extends naturally to higher parameter dimensions, the LP region in Fig. 7 is decomposed 
into triangular regions (blue). An acceptable candidate parameter set (red) is efficiently identifiable by being 
interior to one of the triangles. In contrast, an unacceptable parameter set (black) is not interior to any of the 
triangles and would result in a rollback. 



4.4 Elicitation of Prior 

A valid posterior probability requires multiplication of the likelihood in Eqn. 5 by the prior 
probability of the candidate parameters. In determining an informative prior for / app and 
g syn , the nonlinear nature of the model should not be underestimated. Typically parameters 
interact in highly nonlinear ways to establish the behavior of the system. The different re- 
gions of parameter space can, in principle, be assigned distinct prior probabilities. Another 
approach would be to define a density that varies smoothly and continuously over the pa- 
rameter support. However, in either case the p-dimensional improper integral, taken over 
the support of all parameters, of the joint prior distribution must evaluate to 1.0. If some 
regimes, but not others, are to be considered then an p-dimensional density over a typically 
irregular region must be defined. There is typically little guidance about the shape of such 
regions, in particular such regions may not even be simply connected. 

Our approach regarding elicitation of the prior follows from the typical goals of mathemat- 
ical modeler who wants to know if their model is capable of predicting the response while 
capturing certain key features of interest. Consider that a constant voltage set at the average 
voltage of a neuron firing action potentials may in fact provide a relatively high predictive 
accuracy since it corresponds to the least squares linear regression. However the average 
voltage model, clearly would be unacceptable because it predicts poorly the additional fea- 
tures of amplitude and frequency. We assume researchers will have sufficient expertise to 
know that their data exhibits action-potentials and feel justified in setting to zero the prior 
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probability of any parameters which do not predict such behavior. 

We find that this "NO ROLLBACK" region is a simply connected finite region (having area 
of A w 933.75 /iA ■ mS ■ cm' 4 ). We assume flat uniform prior over this EAS/AAS region. 
This is the appropriate approach when an investigator is interested in inference for a model 
conditioned on a prior belief that the neurons are indeed firing action potentials, but has 
no further prior information about values of J app and g syn . Such flat priors are commonly 
called uninformative priors and are common choice of prior in Bayesian statistical analysis. 
The "no rollback" region of Fig. 7 is assumed to have prior probability equal to 1/A while 
all other parameter values have prior probability zero*. No n- uniform priors over this region 
could be defined which would have the effect to drive parameter estimates towards low/high 
values of J app and g syn . Non-uniform priors are beyond the scope of the present paper but 
are areas for continuing research by our group. 

At each iteration of MCMC, candidate parameters generated by the proposal distribution 
were tested to determine if they resided within the "NO ROLLBACK" region by using the 
MATLAB function pointLocation with the look-up values in Table 2. This approach was 
highly efficient and did not contribute significantly to the computational effort. Candidate 
parameters falling outside the region could be immediately rejected avoiding any further 
computational expense on that iteration. 

4.5 Proposal Distribution 

Metropolis-Hastings sampling in MCMC is distinguished by allowing for the next candidate 
parameters in the chain to be sampled from a proposal distribution. This is convenient 
in applications such as ours where an analytic form for the posterior is either unknown, 
intractable, or computationally inefficient. Typically a proposal distribution is selected which 
i) shares the same support as the parameter(s) of interest and ii) is easy to evaluate. Typically 
so-called location-scale distributions are utilized so that the next guess is drawn from a 
distribution may be centered on the most recent element in the chain. The scale of the 
proposal is often called a mixing value as it controls, in part, how vigorously parameter 
space is traversed. The mixing values may also be adaptively tuned during MCMC. This 
so-called adaptive MCMC is especially useful or even required in nonlinear problems where 
characteristically the sensitivity of the likelihood and/or prior to small changes in parameter 
values can change drastically as the parameter space is traversed. 

Both J app and g syn are physiological parameters. In our application we were only interested in 
excitatory input and positive synaptic conductances. To achieve unconstrained sampling of 
the positive semi-infinite support of these parameters we estimated the log of the parameters' 
values. The proposal distribution for both log-parameter was assumed to be a Gaussian. 
Mixing values (standard deviations) for these proposals was initially assigned to 0.01, but 
then adaptively varied by the operator down to 0.001 after burn-in was determined to have 
occurred. 

*Any positive constant (e.g. 1.0) would serve equally well. 



17 



4.6 MCMC Without Rejection Sampling 



6 
5 
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-2 
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Figure 9. MCMC trace plot 

The blue line is the natural log of the current estimation of the -f app . The green line is the natural log of the 
current estimation of the g S yn- The burn-in phase of this particular trace, where the current estimate takes larger 
steps to find a more likely estimate, ends after 450 iterations. The post-burn-in phase is after this point, where 
the program focuses in on a small ranges of highly likely estimates. 



For purposes of comparison, we now illustrate the MCMC estimation of / app and g syn without 
rejection sampling by ignoring completely the feasibility region illustrated in Fig. 7. In Fig. 9, 
the parameter trace of the MCMC program is shown for one thousand iterations, where the 
moving blue line is the natural log of the current I app estimate and the moving green line is 
the natural log of the current g syn estimate. The flat sections where the guess is unchanged 
for a number of time steps represent a string of rejects where the previous guess, 9 , was 
kept as the most probable estimate. The regions where the trace jumps to a new estimate 
is where a new set of estimates was accepted and replaced the previous. 

The trace can be split into two sections. First is the burn-in phase, where the current 
estimate is jumping relatively far from one set of parameter guesses to another as it comes 
closer to the most likely set of values. The second part is the post-burn-in phase, when the 
estimate limits itself to a very small range and remains more or less stationary around the 
mostly likely sets of values. For Fig. 9, the burn-in phase ends at about 450 time steps. To 
find the best estimate, we take the average of the post-burn- in phase after this point. We 
ran the MCMC routine for five trials in the process outlined above. For known parameters 
and the listed initial guesses, the program finds the data in Table 3. 




iapp 

Bsyn 

psealel 

pscale2 

pleakl 

— pleak2 
mstdl 
mstd2 
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a) Data with I app = 120, Initial guess: 220 b) Data with g syn = 7.5, Initial guess: 1 



Irial 


TV /¥ 

Mean 


CI] J "1 ~v 

btd. Dev. 


1 


238.30 


1.00 


2 


242.53 


1.00 


3 


226.71 


1.00 


4 


240.80 


1.00 


5 


238.11 


1.00 


Overall 


237.29 


1.00 


Percent Error 


97.7% 



Irial 


TV /¥ 

Mean 


btd. Dev. 


1 


0.00 


1.83 


2 


1.42 


1.04 


3 


3.77 


1.00 


4 


1.53 


1.03 


5 


0.00 


2.57 


Overall 


1.35 


1.49 


Percent Error 


82.1% 



Table 3. 7 a pp and gsyn estimation results without rejection sampling 

Means and standard deviations were easily determined as sample statistics calculated over the MCMC trace after 
burn-in. 



Obviously, these estimations are less than satisfactory with high percent errors. For both 
parameters, the program fails to approach the true parameters or deviate from the initial 
guesses. In addition to parameter sets producing EAS/AAS, elements of the MCMC chain 
were also occasionally SS parameter sets. This created flat-line cumulative power and mean 
voltage graphs. While this deviation from the oscillating state of the system was often 
rejected by MCMC, it typically lead to variance components (pscale*, mstd*) being over- 
estimated and if the range searched for more probable values was too small, the chain often 
become trapped - predicting only implausible parameter sets. By letting the program run 
without informed parameter constraints , MCMC produces unrealistic parameter estimates 
for different states. Thus, to find the feasible region where the EAS/AAS solutions exists, 
we look to bifurcation analysis. 



4.7 MCMC With Rejection Sampling 

Five more trials of this MCMC code were run in the same manner described earlier, but 
in combination with rejection sampling. The resulting estimates can be found in Table 4. 
These trials are much more accurate than those without a limited parameter space, as can 
be seen from the dramatic decrease in percent error from Table 3. Both the estimations of 
-f app and g syn moved farther from the initial guess value and gravitated toward the actual 
value. These results are especially satisfying because the initial guesses were so far from the 
actual data values. 
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a) Data with I app = 120, Initial guess: 220 



Irial 


TV /¥ 

Mean 


btd. Dev. 


1 


102.08 


1.03 


2 


109.32 


1.10 


3 


98.79 


1.02 


4 


110.67 


1.07 


5 


100.23 


1.02 


Overall 


104.22 


1.05 


Percent Error 


13.6% 



Table 4. / ;w and </ syn using proposal with rejection sar 

Means and standard deviations were easily determined as 
burn-in. Compared to results without rejection sampling 



b) Data with g syn = 7.5, Initial guess: 1 



Irial 


TV /¥ 

Mean 


btd. Dev. 


1 


8.17 


1.03 


2 


7.62 


1.10 


3 


8.43 


1.03 


4 


7.74 


1.05 


5 


8.34 


1.02 


Overall 


8.06 


1.05 


Percent Error 


7.46% 



npling 

sample statistics calculated over the MCMC trace after 
(see Table 3), these estimates are less biased. 




Time Time 
Figure 10. Predicted voltage vs. time 

The blue line is the true voltage for the data given, and the red line is the voltage for the current parameter 
estimation. 



5 Conclusion 

Estimating parameters with MCMC alone is not enough to procure good results for a dy- 
namical system with different behavioral states. One method we find to be effective for 
the coupled ML model is limiting the parameter region being considered in MCMC using 
rejection sampling. This can be accomplished through first knowing the state of a data set, 
then determining parameter boundaries of that state using bifurcation continuation analysis. 
This removes proposed parameter sets of MCMC that are behaviorally far from the data. In 
unrestrained MCMC, these sets make it difficult for estimates to ever be "accepted," because 
the proposed conditioning statistics are far from that of the data. By applying this method 
to the estimation problem, there is a large decrease in percent error of the estimation param- 
eters, as can be seen in Tables 3 and 4. With the primary goal of estimating / app and g syn , 
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we introduced the concomittant parameters pscale^ pscale 2 , pleak 1? pleak 2 , mstd 2 , mstd 2 
which linked to the parmaters of interest through specific choices modeling the lack of fit of 
the predicted voltage and its cumulative power. As is standard practice in Bayesian latent or 
strucutral equation modeling, these extra parameters effectively represent an hypothesized 
latent covariance relation among the parameters of interest and are considered just as inte- 
gral to model design as the choice of the ML equations or prior distribution. This approach 
proves to be both a convenient and effective method for estimating the coupling strength 
and applied current of two reciprocally coupled ML neurons. 



6 Further Research 

These results cover a small fraction of the potential research that can be done with the 
ML system and parameter approximation. The problem can become increasingly more 
complex with the alteration of initial conditions, addition of dynamical noise, or estimation 
of additional parameters. 

Extension of this method to higher number of parameters is an obvious goal. Triangulation of 
multi-dimensional parameter regions is already possible utilizing N-D Delaunay triangulation 
via MATLAB's delaunayn and tsearchn functions. The challenge therefore is not in determin- 
ing membership of a proposed candidate parameter set within an "NO ROLLBACK" region, 
but rather the higher dimensional parameter continuations need to determine the region to 
be triangularized. In particular, the current version of XPPAUT does not easily permit 
continuation beyond 2D, however, the current version of the AUTO library upon which it 
is based is developing reliable higher dimensional continuation. Automating the process of 
LP and Hopf continuation, possibly through scripting, would obviously be advantageous as 
well. 

In real-world instances of neuron coupling, interference from nearby neuron firings is present 
in the form of dynamical noise. This can be accounted for by setting 5 > 0, in the ML 
model. The presence of dynamical noise can cause the system to experience MMO, or the 
stochastic switching between states (shown in Fig. 11). The parameters of this ML variant 
may still be able to be estimated, but alterations in the estimation strategy will be necessary. 
For instance, 5 could be estimated in addition to J ap p and g syn . Suitable adjustments to the 
rejection region could be made to account for stochastic perturbations of I app . 
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Figure 11. Mixed mode oscillation 

In the presence of small noise (5i = 82 = 0.7) and weak coupling (g syn = 0.15) the reciprocally coupled ML 
model with I app = 95 predicts low amplitude oscillations disrupted by large amplitude excursions. 



The parameters of the ML system may also be estimated by conditioning on other statistics. 
In this study we used cumulative power and mean of the voltage to run MCMC, but this 
could be altered. For example, conditioning on the difference of cumulative powers of the two 
voltages was considered, but has not been implemented. This approach may better account 
for synchrony, or lack thereof, between the two neurons than calculating the likelihoods for 
each neuron independently. 

Appendix 

A.l Initial Conditions and Parameters 

Unless indicated otherwise, all tests run on the MCMC system had the following initial 
conditions and parameter values. A biological interpretation of these parameters can be 
found in Table 1. 



Variable 


Initial Value 


Variable 


Initial Value 


Vl 


-20 mV 


v K 


-84 mV 


v 2 


20 mV 


vl 


-60 mV 


Wi 


0.3 


Vsyn 


70 mV 


w 2 


0.5 


Vn 


-1.2 mV 


Si 


0.2 


V22 


18 mV 




0.1 


v 3 


2.0 mV 


C 


20 fiF cm 2 




30 mV 




4.0 mS/cm 2 


T 


8 ms 


9k 


8.0 mS/cm 2 


v B 


5.0 mV 


9l 


2.0 mS/cm 2 


v t 


15 mV 




120 mS/cm 2 





0.04/ms 
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A. 2 Cumulative Power of Periodic Functions is 0(t) 
Let m(t) have Fourier series representation, 



m(t) = ^ + ^2 a k cos(27T0 fe t) + b k sin(27r0 fe t) 



k=i 



then, it it to be shown that 

p(t) = c-t + o(i) 

Proof is By Induction and is adapted from Quinn (2011). In the base case (n = 1), 

9 (t) e 0(1) 
. * . 



P(t) = [87r 4 4 (a 2 + 6?)] t +27r 3 a50?sin(47r0iO 

- 2^b\(f)\ sin(4vr0it) 

+ 167r 4 ai&i0i cos 2 (27r0!t) 

— 167rai6i0 4 



Then it is supposed that for n — 1 we have, 



P{t) = 1^ I 4:ir 2 a k </) 2 k cos(27T0 fc t) + 4ir 2 b k <f) 2 k sin(27r0 fe t) j dt 



K k=l 

and it remains to be shown that 
"t ( n 



n-l 



t + g n -i(t) 



k=i 



P(t) = jT 4vr 2 a fc 0^ cos(27T0 fc t) + 4tt 2 6 A; 0^ sin(2vr0 fe t) j rft 

Next, collecting the n th terms from the summation yields, 

P W=/ y"47r 2 a fc 2 cos(2vr0 fe t) 
Jo \k=i 

+ Air 2 b k (j) 2 k sin(27T0A;t) + 47r 4 a n 2 cos(27T0 n t) 



i + 9n(t) 



k=l 



+47r 4 6 n 2 sin(2vr0 n t) 
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Expanding the square in the previous result gives, 



P(t) = ^2 cos(27r0 fc t) + A-K 2 b k (pl sin(27T0 fe t) j dt 



t I n 



+ 2 j ^ a k4>l cos(27T0 fc t) + Aix 2 b k 4>l sin(27T0 fc t) 

■ (47r 2 a„0 2 cos(27T0„t) + 47r 2 6 n (/) 2 sin(27T0 n t)) (it 
+ / (47r 2 a n 2 t cos(2?r0 n t) + 4tt 2 6„0 2 sin(27T0 n t)) 2 di 



The first term is recognized as the induction hypothesis and so has the form C ■ t + g n -i(t). 
The integrals of the remaining terms are evaluated with the extensive use of trigonometric 
identities. 



P{t) = 87T 4 J2 <t>U°-k+K) 
k = l 

n — 1 

+ ^ a k a n4>\^n 

k = l 

n—1 

- 8TT 3 b k a n rf>^4 > ^ l 
k = l 
n — 1 

+ 23 S7 ' 3b k a n<t>k<t>' n 



fc = l 
n — 1 

S7r3a kbn4>l<l>l, 

k = l 
n — 1 

^2 Sn 3 b k bn<l>l4>n 

k=l 



t + Sn-lM 



"sin(27r(0 fc 


-*„)«) + 


sin(2ir(0 fc 


+ «„)*)" 


<*fc - 




<^fc + 




cos(2ir(0 fc 


- 4>n)t) 


cos(2ir(0 fc 


+ <*>„)*)' 


4>k - 




4>fc + 


tfn 


1 


1 






4 

. <t>k — 4>n 


4>k + 4>n _ 






cos(2ir(0 fc 


-4>n)t) ! 




+ <*>„)*)' 


4>k - 




4>fc + 




1 


1 






_<t>k — <t>n 
sin(2?r(0 fc 


4>k + 4>n _ 
-4>n)t) 


sin(2ir(0 fc 




0fc - 


<t>n 


0fc + 





+8ir 4 <ji*bjt - 2tt 3 6^* sin(47T0„i) 



The bold terms may be combined with the leading term of the induction hypothesis raising 
the upper bound of the summation from n — 1 to n. The remaining terms, only containing 
t as arguments of sines and cosines, can be merged with g n -i{t) from the induction 
hypothesis. Calling this merger g n (t) completes the induction. Cumulative power has been 
written in the desired form 



P(t) = 



at + bl) 



t + g n (t) 



= C-t + 0{l) 



Since clearly C -t e 0{t) and g n (t) G 0(t) it follows that their sum P(t) £ 0(t) proving 
the desired result. 
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